#ifndef _definition_h_
#include "lib.h"
#define _definition_h_
#endif
//contructor
double** contructor(int row, int col){
    double** result = new double* [row];
    for(int e=0;e<row;e++){
        result[e] = new double[col];
    }
    for(int i=0;i<row;i++){
        for(int j=0;j<col;j++){
            result[i][j]=0;
        }
    }
    return result;
}
//**************************************************
//Get_xyz_data

double** readFile( int* num, char* file) {

    using namespace std;
	ifstream fin(file);
    string line;
	if (!fin) {
        cout <<"Error reading " << file << endl;
		return NULL;

	};
//	*num = 1;
//    for (int i=0;i<R;i++){
//
//        for (int j=0;j<100;j++){
//            fin >> data[i][j];
//            *num = *num+1;
//            if (fin.eof()){
//                    *num = *num/R;
//                    fin.close();
//                    return 0 ;
//            }
//        }
//    };
    double value;
    int i=0,j;
    while (getline(fin,line)) {
            *num=j;
            std::stringstream iss(line);
            j=0;
       while(iss >> value)
    {
        j++;
    }
    i++;
    }
    i=0;
    double** data=new double*[R];
    for(int e=0;e<R;e++){
        data[e] = new double[*num];
    }
    fin.close();
    ifstream f_in(file);
    while (getline(f_in,line)) {
        std::stringstream iss(line);
        j=0;
       while(iss >> value)
    {

        data[i][j]=(value);
        j++;
    }
    i++;
    }

	f_in.close();
	return data;
}
//count element
//int count_element(double** data){
//      int i=0,j=0;
//     while(data) {
//           //cout << data[i][0] << endl;
//        data++;
//        i++;
//     }
//
//     return i;
//     }
//Find max in data
double** max_data(double** data,int num){
    double **max_dt;
    max_dt = contructor(1,num);
    for(int i=0;i<num;i++){
        for(int j=0;j<R;j++){
            if (data[j][i]>max_dt[0][i]) max_dt[0][i]=data[j][i] ;
        }
    }
    return max_dt;

}

//find min in data
double** min_data(double** data,int num){
    double **min_dt;
    min_dt = contructor(1,num);
        for(int j=0;j<num;j++){
            min_dt[0][j]=10000 ;
        }
    for(int i=0;i<num;i++){
        for(int j=0;j<R;j++){
            if (data[j][i]<min_dt[0][i]) min_dt[0][i]=data[j][i] ;
        }
    }
    return min_dt;

}
//sum x,y,z data in a matrix
double sum_matrix(double** data,int ro,int co){
    double result=0;
    for (int i=0;i<ro;i++)
        for(int j=0;j<co;j++){
            result=result+data[i][j];
        }
    return result;
}
//Sub matrix
double** sub_matrix(double** data1, double** data2,int ro,int co){
     double** result;

     result = contructor(ro,co);
    for (int i=0;i<ro;i++)
        for(int j=0;j<co;j++){
            result[i][j]=data1[i][j]-data2[i][j];
        }
    return result;
}
//Mul matrix
double** mul_matrix(double** data1, double** data2,int num){
    double** result;
    result = contructor(R,num);
    for(int i=0; i<R;i++){
        for(int j=0;j<num;j++){
            result[i][j]=0;
            for(int x=0; x<R;x++){
                result[i][j]= result[i][j] + data1[i][x]*data2[x][j];
            }
        }
    }
    return result;
}
//****************************
//get_point_centroids >OK
double** get_point_centroids(data_t data,int K, int D){
    double** mean;
    mean = contructor(R,D);
    for (int n=0;n<R;n++){
        for (int i=0;i<data.num;i++ ){
            for (int j=0;j<D;j++){
                if(j==0)
                    mean[n][j] = mean [n][j]+data.data_x[n][i];
                else if (j==1)
                    mean[n][j] = mean [n][j]+data.data_y[n][i];
                else if (j==2)
                    mean[n][j] = mean [n][j]+data.data_z[n][i];

            }

        }
        for (int j=0;j<D;j++){
            mean[n][j]=mean[n][j]/data.num;
        }
    }
    double** centroids;
    centroids= kmeans(mean,K);
    return centroids;
}
//***************************************
//kmeans > OK 0.4
//nbCluster = N :states
double**  kmeans(double** data, int nbCluster){
    int data_dim, nbdata;
    double** data_min,**data_max,**data_diff;
    data_dim = D;
    nbdata = R;
    double** centroid,**diff_Position;
    centroid = contructor(nbCluster,data_dim);
//init the centroids randomly
    data_min = min_data(data,D);
    data_max = max_data(data,D);
    data_diff = sub_matrix(data_max,data_min,1,D);
//every row is a centroid
   srand((unsigned)time(NULL));
//data_diff
    for(int i=0;i<nbCluster;i++){
        for (int j=0;j<data_dim;j++){
            centroid[i][j]=0.4;

        }
    }
    for(int i=0;i<nbCluster;i++){
        for (int j=0;j<data_dim;j++){
            centroid[i][j] = centroid[i][j]*data_diff[0][j];
            centroid[i][j] = centroid[i][j]+data_min[0][j];
        }
    };
//end init centroids
//yes>>
    double pos_diff = 1,min_diffs,curAssignment,diff2cs;
    double min_diff[3],diff2c[3];
    int *assignment=new int();
    double pointsInCluster[nbCluster];
    //main loop
    while (pos_diff > 0.0){
        delete assignment;
        int cur=0;
        int *assignment=new int();

//E-Step
    for (int d=0;d<R;d++){
        for(int i=0;i<D;i++){
            min_diff[i]=data[d][i] - centroid[0][i];
        }
        min_diffs=0,diff2cs=0;
        for(int i=0;i<D;i++){
            min_diffs = min_diffs + min_diff[i]*min_diff[i];
        }
        curAssignment = 1;
        for (int c=1;c<nbCluster;c++){
            for(int i=0;i<D;i++){
                diff2c[i]=data[d][i] - centroid[c][i];
            }
            diff2cs=0;
            for(int i=0;i<D;i++){
            diff2cs = diff2cs + diff2c[i]*diff2c[i];
            }
            if(min_diffs>=diff2cs){
                curAssignment = c;
                min_diffs = diff2cs;
            }
        }
        //assign the d-th dataPoint
        assignment[cur]=curAssignment;

        cur++;

    }

    //for the stoppingCriterion
    //oldPositions = centroid;
    double** oldPositions;
    oldPositions = contructor(nbCluster,data_dim);
    for(int i=0;i<nbCluster;i++){
        for (int j=0;j<data_dim;j++){
            oldPositions[i][j] = centroid[i][j];
        }
    };
    //M-Step
    //recalculate the positions of the centroids
        for(int i=0;i<nbCluster;i++){
        for (int j=0;j<data_dim;j++){
            centroid[i][j] = 0;
        }
    };
    for(int i=0;i<nbCluster;i++){
        pointsInCluster[i]=0;
    }
    for (int d=0;d<cur;d++){
        for(int i=0;i<data_dim;i++){
            centroid[(assignment[d])][i]=centroid[(assignment[d])][i]+data[d][i];
        }
        pointsInCluster[(assignment[d])]=pointsInCluster[(assignment[d])]+1;
    }


    for (int c=0;c<nbCluster;c++){
        if(pointsInCluster[c] != 0){
            for(int i=0;i<data_dim;i++){
                centroid[c][i]=centroid[c][i] / pointsInCluster[c];
            }
        }
        else{
                // set cluster randomly to new position
            for(int i=0;i<data_dim;i++){
                //centroid[c][i]=(((double) (rand() % 10000)/10000) * (data_diff[0][i]))+data_min[0][i];
                centroid[c][i]=(0.4* (data_diff[0][i]))+data_min[0][i];
            }
        }
        //stoppingCriterion

    }

    diff_Position=sub_matrix(centroid,oldPositions,nbCluster,data_dim);
    for (int i=0;i<nbCluster;i++){
        for(int j=0;j<data_dim;j++){
            diff_Position[i][j]=diff_Position[i][j]*diff_Position[i][j];
        }
    }
    pos_diff = sum_matrix(diff_Position,nbCluster,data_dim);

}

 //   *data_diff[0]=
//    for(int i=0;i<D;i++){
//        cout<<data_diff[0][i]<<'a';
//    }

     return centroid;

}

//get_point_clusters >OK
double** get_point_clusters(data_t data, double** centroids,int D){
    double** XClustered;
    double temp[N];
    double temp1,temp2,temp3;
    XClustered = contructor(data.num,R);
    for(int n=0;n<R;n++){
        for(int i=0;i<data.num;i++){
            for(int e=0;e<N;e++){
                temp[e]=0;
            }
            for(int j=0;j<N;j++){
                if(D==3){
                    temp1=(centroids[j][0] - data.data_x[n][i])*(centroids[j][0] - data.data_x[n][i]);
                    temp2=(centroids[j][1] - data.data_y[n][i])*(centroids[j][1] - data.data_y[n][i]);
                    temp3=(centroids[j][2] - data.data_z[n][i])*(centroids[j][2] - data.data_z[n][i]);
                    temp[j]= sqrt(temp1+temp2+temp3);
                }
            }
            double min_temp=1000000;
            int I_min=0;
            for (int j=0;j<N;j++){
                if(min_temp>temp[j]){
                    min_temp=temp[j];
                    I_min=j;
                }
            }
            XClustered[i][n]=I_min;
        }
    }
    return XClustered;
}
//prior_transition_matrix >OK
double** prior_transition_matrix(int K,int LR){
    //P=(1/LR)*eye(K)
    double** P;
    P=contructor(K,K);
    for(int i=0;i<K;i++)
        for(int j=0;j<K;j++){
            if(i==j){
                P[i][j]= (double) 1/LR;
            }
        }
    for(int i=0;i<(K-(LR-1));i++){
        for(int j=1;j<=(LR-1);j++){
            P[i][i+j] = (double) 1.0/LR;
        }
    }
    for(int i=(K-(LR-2)-1);i<K;i++){
        for(int j=1;j<=(K-i);j++){
            P[i][i+(j-1)] = 1.0/(K-i);
        }
    }
    return P;
}
//*************************************************************************************************
//dhmm_numeric
double** dhmm_numeric(double** X, double** P, int *bins, int K, int cyc, double tol){
//% function [E,P,Pi,LL]=dhmm_numeric(X,alphabet,K,cyc,tol);
//%
//% simple Hidden Markov Model - variable lengths and discrete symbols
//%
//% X - N x p array of Xs, p sequences
//% bins - possible bins of Xs
//% K - number of states (default 2)
//% cyc - maximum number of cycles of Baum-Welch (default 100)
//% tol - termination tolerance (prop change in likelihood) (default 0.0001)
//%
//% E - observation emission probabilities
//% P - state transition probabilities
//% Pi - initial state prior probabilities
//% LL - log likelihood curve
//%
//% Iterates until a proportional change < tol in the log likelihood
//% or cyc steps of Baum-Welch
    int _D=1;
    int num_bins = N;
    double epsi = 1/10000000000;
    // N = size (X,1);
    double *T = new double();
    for (int n = 0;n<N;n++){
        T[n]= R;
    }
    //TMAX=max(T);
    int TMAX = R;
    //cout <<"*****************************************************************************";
    //cout<<"Training
    //E = (0.1*rand(num_bins,K)+ones(num_bins,K))/num_bins;
    double** E;
    E = contructor(num_bins,K);
    for(int i=0;i<num_bins;i++){
        for(int j=0;j<K;j++){
            E[i][j] = (0.1*0.4 +1.0)/num_bins;
        }
    }
    double* sum_E;
    for(int i=0; i<K;i++){
        sum_E[i] = 0;
    }
    for(int i=0; i<K;i++){
        for(int j=0;j<num_bins;j++){
            sum_E[i]=sum_E[i]+E[j][i];
        }
    }
    double** Z ;
    Z= contructor(num_bins,K);
    for(int i=0; i < 8 ; i++){
//        for(int j=0 ; j<K ; j++){
//            //Z[0][0] = 0;
//        }
    }
//    for(int i=0;i<num_bins;i++){
//        for(int j=0;j<K;j++){
//            Z[i][j]=E[i][j] / sum_E[j];
//        }
//
//    }









    return E;
}
